Critical behavior of the Widom-Rowlinson mixture: coexistence diameter and order 

parameter 
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The critical behavior of the Widom-Rowlinson mixture [J. Chem. Phys. 52, 1670 (1970)] is studied 
in d = 3 dimensions by means of grand canonical Monte Carlo simulations. The finite size scaling 
approach of Kim, Fisher, and Luijten [Phys. Rev. Lett. 91, 065701 (2003)] is used to extract the 
order parameter and the coexistence diameter. It is demonstrated that the critical behavior of the 
diameter is dominated by a singular term proportional to i 1_a! , with t the relative distance from 
the critical point, and a the critical exponent of the specific heat. No sign of a term proportional to 
t 2/3 could be detected, with /3 the critical exponent of the order parameter, indicating that pressure- 
mixing in this model is small. The critical density is measured to be pa 3 — 0.7486 ± 0.0002, with a 
the particle diameter. The critical exponents a and /3, as well as the correlation length exponent v, 
are also measured and shown to comply with d = 3 Ising criticality. 

PACS numbers: 02.70.-c, 05.70.Jk, 64.70.Fx, 64.60.Fr 



I. INTRODUCTION 

The Widom-Rowlinson (WR) mixture^ is a simple 
model of a fluid exhibiting phase separation. The model 
consists of A and B particles that interact via simple pair 
potentials: the AA and BB pair interaction is ideal, while 
AB pairs interact via a hard-core potential of diameter 
a. Upon increasing density, the WR mixture phase sep- 
arates into an A-rich and B-rich phase. In d = 3 dimen- 
sions, computer simulations agree that the corresponding 
universality class is that of the three-dimensional (3D) 
Ising mo defies 5 .. There is, however, some variation in 
the reported estimates of the critical density. 

In order to probe fluid criticality using computer sim- 
ulation, high quality data are required. The latter are 
typically generated using Monte Carlo (MC) methods, 
and considerable effort has been devoted to develop effi- 
cient MC schemes. In Ref. 0, for example, a MC cluster 
move is described for the WR mixture that is (nearly) 
free of critical slowing down. However, as pointed out in 
Ref. H, this type of move cannot be used to obtain the 
coexistence curve, which may be problematic if one is in- 
terested in measuring, say, the critical exponent of the 
order parameter. In Ref. H, therefore, a different clus- 
ter move is formulated, based on Ref. |fj, which not only 
gives access to the coexistence curve, but is also rejec- 
tion free. Recently, the latter approach was generalized 
to continuous potentials^. 

However, in addition to efficient MC sampling, of at 
least equal importance (if not more) is the finite size scal- 
ing (FSS) algorithm used to extrapolate the simulation 
data to the thermodynamic limit. FSS is essential be- 
cause the correlation length diverges at the critical point, 
and thus the true thermodynamic limit is never captured 
in a finite simulation box, no matter how efficiently it is 
simulated. For fluids, recently proposed unbiased FSS 
algorithms formulated in the grand canonical ensemble 
seem particularly powerfu l 8 ^ 10 ' 11 . The latter algorithms 



are unbiased in the sense that no prior knowledge of the 
universality class is required: the critical point of the 
transition, as well as some of the critical exponents, are 
an output. These unbiased algorithms were used, for ex- 
ample, to resolve the universality class of the hard-core 
square-well (HCSW) fluid and the restricted primitive 
electrolyte, both of which were shown to exhibit 3D Ising 
critical behavios2ii£. 

Unfortunately, it is not obvious how the MC cluster 
moves for the WR mixture generalize to the grand canon- 
ical ensemble. Grand canonical cluster moves for mix- 
tures seem less common, but some have been presented 
in the literature 13 ' 14 . Of these, the MC move of Ref. El 
a generalization of its canonical variant 15 , is readily ap- 
plicable to the WR mixture. By using the MC move of 
Ref. Q, a FSS analysis of the WR mixture using the 
above mentioned unbiased algorithms thus becomes pos- 
sible. This, consequently, is the aim of the present work. 
Of particular interest is the coexistence diameter, whose 
critical behavior is governed by a very weak singularity 
that is challenging to extract from simulation data. Note 
that, at the time of writing, the approach of Ref . seems 
to be the only FSS algorithm available to extract the 
coexistence diameter correctly from simulation data. A 
correct description of the latter is required in order to reli- 
ably estimate the critical density^. Since the coexistence 
diameter of the WR model has not received much atten- 
tion in previous simulations, the present grand canonical 
approach is certainly warranted. 



II. SIMULATION METHOD 

In the grand canonical ensemble, the volume V and 
the fugacity z\ (zb) of species A (B) are fixed, while 
the particle numbers Na and N-q fluctuate. The thermal 
wavelength is set to unity such that the fugacity z a di- 
rectly reflects the number density N a /V a pure phase of 
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a particles would have (recall that such a phase is simply 
an ideal gas). In what follows, all remaining length scales 
are expressed in terms of the hard-core diameter a. The 
crucial quantity is the finite-size grand canonical distri- 
bution Pl(Na, -/Vb|za, zb), defined as the probability of 
observing a system containing Na particles of species A 
and ,/Vb particles of species B, at fugacities za and Zb, 
with L the lateral dimension of the cubic simulation box 
(the use of periodic boundary conditions is assumed) . 

The distribution is obtained numerically in a grand 
canonical MC simulation via the insertion and removal 
of particles. To simulate efficiently, the cluster move of 
Ref. 1 1 ll is used; Pl{Na, Nb\za, zb) is then obtained sim- 
ply by maintaining a histogram. To overcome the free 
energy barrier separating the phases, a biased sampling 
scheme is also implemented^. Here, the simulation is di- 
vided into distinct intervals (called windows) each span- 
ning a single A particle, while Nb is allowed to fluctuate 
freely inside each window (the choice for A or B is ar- 
bitrary). The windows are then sampled separately and 
successively. CPU time is divided such that the number 
of "sweeps" per window is the same for all windows. In 
this work, we say that a sweep has passed when a given 
population of particles has completely been replaced or 
updated by new ones. This is in contrast to the more 
common approach of keeping the number of attempted 
MC moves per window fixed. The latter approach, how- 
ever, is less appropriate for grand canonical simulations 
since the acceptance rate is typically density dependent. 
Per window, approximately 1800 sweeps are generated. 
To obtain a single distribution, an investment of around 
7 CPU hours for a small system (L = 8), and 270 hours 
for a large system (L = 13) is required. In order to per- 
form the subsequent FSS analysis, Pl{Na, Nb\za, zb) is 
measured for system sizes L = 8 — 13 at fugacities ranging 
from close to the critical point to well into the coexistence 
region. Estimates of properties at intermediate fugacities 
are obtained using the multiple histogram method 1 ^. 



III. RESULTS 

As mentioned before, the WR mixture exhibits phase 
separation into an A-rich and B-rich phase. Note that 
the B-rich phase may equally well be regarded as being 
poor in A species. In this sense, then, phase separa- 
tion is analogous to liquid-vapor coexistence: the A-rich 
phase being the liquid, the A-poor phase being the va- 
por, and the fugacity of the B particles being inverse 
temperature (again, the choice for A or B is arbitrary). 
A natural definition of the order parameter is therefore 
A = (pl — Pv)/2, with pl the number density of A parti- 
cles in the A-rich phase, and py the number density of A 
particles in the A-poor phase. Close to the critical point, 
the order parameter is expected to scale as A cx t 13 , with 
t = zb/ zb,ct — 1 the distance from the critical point, and 
zb.ci- the critical "inverse temperature" . Similarly, the 
coexistence diameter can be written as D = (pl + Pv)/2- 



The critical behavior of the latter is given bjsi^ 

D = pa.c, (1 + A 2l3 t 213 + A 1 _ a t 1 ~ a + A ± t) , (1) 

with pA,cr the number density of A particles at the critical 
point, and non-universal amplitudes Ai. For the 3D Ising 
universality class, appropriate exponent values are (3 « 
0.326 and a w 0.1092&. 

A. Order parameter 

To extract the order parameter, the FSS algorithm of 
Ref. |9j is used (for a more detailed description of the al- 
gorithm, Ref. [l0| is also highly recommended) . The al- 
gorithm requires as input the grand canonical distribu- 
tion Pl(Na, Nb\za, zb) for at least three system sizes 
L. Here, five system sizes L = 9,10,11,12,13 are in 
fact used. Starting with zb significantly above its crit- 
ical value, the cumulant ratio (m 2 ) /(m 4 ) is plotted as 
function of the average number density (Na) /V, with 
m = Aa — (Na) (note that this plot is parameterized by 
the fugacity of the other species za)- The resulting curve 
will reveal two minima, located at p~ and p + , with re- 
spective values Q~ and Q + at the minima. Defining the 
quantities Q min = (Q + + <2~)/2, x = Q mhl ln(4/eQ mi „), 
and y = (p + — p~)/(2A), the points (x,y) from the dif- 
ferent system sizes should, in the limit far away from the 
critical point, collapse onto the line y = 1 + x/2. Re- 
call that A is the order parameter in the thermodynamic 
limit at the considered fugacity zb, precisely the quan- 
tity of interest, which may thus be obtained by fitting 
until the best collapse onto 1 + x/2 occurs. In the next 
step, zb is chosen closer to the critical point, the points 
(x,y) are calculated as before, but this time A is cho- 
sen such that the new data set joins smoothly with the 
previous one, yielding an estimate of the order parame- 
ter at the new fugacity. This procedure is repeated all 
the way to the critical point, where A vanishes, leading 
to an estimate of the critical fugacity zb.ct- Moreover, 
the procedure also yields y as function of x. The latter 
scaling function is universal within a universality class, 
and for the HCSW fluid can be found in Ref. |9j Since 
the WR mixture belongs to the same universality class, 
a similar curve should be found. The latter is verified in 
Fig. 2] which shows y as function of x obtained in this 
work, compared to the result of Ref. The agreement 
is very reasonable. From the vanishing of the scaling 
function, at x c = 0.280, an unbiased estimate of the crit- 
ical fugacity z B , C r = 0.93791 ± 0.00004 is obtained. Note 
that x c is universal within a universality class. The esti- 
mate reported here compares favorably to x c = 0.286 ob- 
tained for the HCSW fluid^, and x c = 0.296 obtained for 
the 3D Ising modelSi, providing additional confirmation 
that these systems belong to the same universality class. 
Shown in Fig. [21 on double logarithmic scales, is the or- 
der parameter A of the WR mixture as function of the 
distance from the critical point t, where the above quoted 
estimate of zb.ct was used. The resolution of the present 
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data is such that A can be resolved down to t RJ 5 x 10~ 5 . 
By fitting the lowest few points in Fig. [3 to A oc t@ , the 
critical exponent is measured to be (3 ~ 0.322 ± 0.008, 
which is certainly compatible with the accepted 3D Ising 
value. 

B. Coexistence diameter 

To extract the coexistence diameter, the FSS algorithm 
of Ref. fll] is used. The algorithm is similar in spirit to 
the previous one, in the sense that it generates a scaling 
function y — f(x), starting with data obtained well away 
from the critical point, and then recursively working its 
way down toward criticality. In Fig. the scaling func- 
tion of the diameter for the WR mixture thus obtained 
is shown, where, as before, five system sizes L = 9 — 13 
were used. For x — ► 0, this function is expected to ap- 
proach y — x/2, which indeed it does. In contrast to the 
order parameter, however, the scaling function of the co- 
existence diameter is not universal— Therefore, a direct 
comparison to scaling functions of other systems cannot, 
in general, be carried out. Nevertheless, for systems with 
negligible pressure mixing, such as the HCSW fluid and 
presumably also the WR mixture, the scaling function is 
expected to be well described by the approximant 11 

e,(x) = d [l-(l-x) 1 -" (2) 
1 + SiSB + s 2 x 2 + S 3 X 3 
1 + ttx + t 2 x 2 + t 3 x 3 J ' 

with x = x/x c , t\ = s% — 1 + a + x c /2Ci, and critical 
exponent a rj 0.109. A fit to the WR data of Fig. EH 
shows that this is indeed the case, with explicit parameter 
values Q = 0.429, x c = 0.175, s x = 4.50, s 2 = -5.72, 
,s 3 = 0.12, ti = 3.81, t 2 = -9.08 and t 3 = 4.25. These 
values are remarkably consistent with estimates quoted 
in Ref. O for the HCSW fluid. Note that ei(x) becomes 
singular close to x c , implied by the (1 — x) l ~ a factor in 
Eq.JSJ. The latter would yield a vertical tangent in the 
plot, at the arrow in Fig. EI The present simulation data, 
however, seem not to extend close enough to the critical 
point to reach this regime. 

The critical behavior of the coexistence diameter D is 
shown in Fig. 0] where ze.cr = 0.93791 obtained in the 
previous paragraph was used. In order to facilitate the 
comparison to other work, 2 x D is actually plotted. Sym- 
metry considerations ensure equal numbers of A and B 
particles at criticality, such that the overall critical num- 
ber density equals p CI = 2 j OA,cr, which is the quantity 
usually quoted in the literature. A fit to the asymptotic 
expansion of Eq.QJ yields p CI = 0.7486 ± 0.0002, where 
the error reflects the variation stemming from the range 
over which the fit is performed (repeating the entire anal- 
ysis leaving out the smallest system size yields a similar 
result). The corresponding amplitudes read as A 2 p rj 0, 
A^ a = 2.76±0.07 and A 1 rj -1.27±0.09, implying that 
the singular behavior is dominated by i 1- ". This, in com- 
bination with the observation that the scaling function is 



well described by ei(x), confirms that pressure mixing 
in the WR mixture is small. The curvature of the di- 
ameter close to the critical point thus reflects the t 1_a 
singularity. Since 1 — a is close to unity, and the magni- 
tudes of Ai_ Q and A\ are similar, the curvature is hard 
to see in Fig. 0] The singular behavior of the diameter 
can be visualized nevertheless by plotting the "inverse 
temperature" derivative k = 2dD/dt instead. In case of 
singular behavior, k is expected to diverge when t — > 0, 
see Eq.Q. Though not very precise, this procedure even 
allows for an unbiased measurement of the exponent a. 
The result is summarized in the inset of Fig. 01 which 
shows k as function of t, where again zb. ct = 0.93791 in 
t was used. The divergence is clearly visible. By fitting 
the data to the form k = a\t~ a +a 2 , with fit parameters 
a, and a, the specific heat exponent is measured to be 
a = 0.11 ± 0.02, which is surprisingly close to the 3D 
Ising value. 



C. Cumulant intersections 



For the sake of completeness, and also to check the 
consistency of the results obtained so far, the critical 
fugacity zb, ct is measured again, but this time around 
using the cumulant intersection approach 22 . As was 
shown by Binder—, the (for example) first order cumu- 
lant Ui = (rn 2 ) / (\m\) 2 becomes system-size independent 
at the critical point. Plots of Ui as function of zb for 
different system sizes are thus expected to show a com- 
mon intersection point, leading to an unbiased estimate 
of the critical fugacity. Moreover, the cumulant value 
Q c at the intersection point is universal, dependent only 
on the universality class. Shown in Fig. [S] is the result 
of this procedure, on a rather fine scale. The resulting 
estimate reads as zb.ct = 0.9379 ± 0.0004, where the er- 
ror reflects the scatter in the intersection points. The 
latter is fully consistent with the previous, more pre- 
cise value, 2 B , cr = 0.93791 ± 0.00004 (arrow in Fig.0). 
For the critical value of the cumulant Q c ss 1.223 is ob- 
tained. This value compares quite favorably to the esti- 
mate Q c = 1.2391 ± 0.0014 obtained in large-scale simu- 
lations of the 3D Ising lattice models, deviating from it 
by less than 2%. The inset of Fig.|S]shows the slope of the 
cumulant Y\ = |dJ7i/dzs| at the critical fugacity, as func- 
tion of the system size L. It is expected that Y% oc L x l v 
with v the critical exponent of the correlation length. Al- 
though there is some scatter in the intersection points, 
the cumulant slopes seem rather constant over the range 
of Fig. |31 and so it is expected that v can be obtained 
quite reliably nevertheless. Indeed, by performing a fit to 
the data in the inset of Fig. [31 the exponent is measured 
to be v w 0.630 ± 0.005, in excellent agreement with the 
accepted 3D Ising value v\ s rj 0.630— 
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IV. DISCUSSION AND SUMMARY 

In this work, the critical behavior of the WR mixture 
was investigated using grand canonical MC simulations 
and unbiased FSS algorithms. As expected, the univer- 
sality class of the transition is that of the 3D Ising model. 
This was demonstrated by direct measurements of the 
exponents a, (3 and v. Substantial indirect evidence has 
also been provided, by comparing the scaling function 
of the order parameter and the coexistence diameter to 
those of the HCSW fluid, as well as via the universality 
of Q c at the cumulant intersection point. 

The critical density obtained in this work 0.7486 ± 
0.0002 can be compared to other simulations. Shew 
and Yethiraj report 0.762 ± 0.016 2 using semigrand 
simulations^^, while Johnson et al. (JEA) obtained 
0.748 ± 0.002 3 . More recent estimates are due to Gozdz, 
0.759 ± O.OISH using the Bruce- Wilding field mixing 
technique 2 ^, and Buhot 0.7470 ± 0.0008 5 . Of the above, 
JEA and Buhot are very close to the value reported here, 
see Fig. The estimate of Gozdz, however, is higher. A 
possible explanation is the adapted FSS algorithm in the 
latter, which seems to overestimate the critical density in 
some cases^. JEA also report an estimate of the critical 
fugacity ze.cr ~ 0.9403 for the largest system size con- 
sidered by themi, but without a systematic FSS analysis 
of this quantity. This overestimates the present value 
significantly. Interestingly, these authors observe an in- 
crease of zs.cr with system size, in disagreement with the 
present work. 

Buhot, by using rejection-free cluster MC moves, is 
able to simulate impressively large systems^, up to L = 
100, which exceeds the typical system size of the present 
investigation by about one order of magnitude. The im- 
proved accuracy of the critical density obtained in this 
work may therefore seem surprising. It should be empha- 
sized, however, that critical phenomena are most conve- 
niently studied in terms of a field variable, such as tem- 
perature or, in the case of the WR mixture, the fugacity. 
The ensemble used by Buhot, as well as the semigrand 



ensemble, do not have access to the fugacity. Instead, in 
these ensembles, the critical point is approached by vary- 
ing the overall density p = (Na + Nb) /V. This somewhat 
restricts the investigation of critical phenomena because 
for every considered p, an explicit simulation needs to 
be carried out. In grand canonical simulations, on the 
other hand, one has access to the particle fugacities. This 
facilitates the extrapolation of simulation data obtained 
at one set of fugacities to different values via histogram 
reweighting 18 . Clearly, the investigation of subtle effects, 
such as the critical behavior of the coexistence diameter, 
is not really feasible without such extrapolation methods. 
Note that the behavior depicted in Fig.^is not simply an 
"artifact" of the grand canonical ensemble. In the semi- 
grand ensemble, for example, the singular behavior of the 
diameter leads to a renormalization of the critical expo- 
nent /3^i. Assuming negligible in Eq. Q , one obtains 
p/Pcr — 1 oc t x ~ a close to the critical point. Combining 
this with the critical power law of the order parameter 
A oc t@ and eliminating t, yields A oc (p/p a - ~ l)' 3 , with 
renormalized exponent (3* — (3/(1 — a). The latter renor- 
malized exponent has been confirmed experimentally 28 , 
and should, in principle, also show up in the WR mix- 
ture when, as mentioned above, the critical point is ap- 
proached by varying p. 

Needless to say, the WR mixture has also been stud- 
ied by theoretical means, using for example density func- 
tional theorjiSS, and integral equations 2 *^. These inves- 
tigations, however, all pertain to the mean-field level. As 
such, quantitative agreement with computer simulations 
close to criticality is not to be expected, and a compari- 
son is consequently not carried out. 
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V. FIGURE CAPTIONS 

FIG. 1: Scaling function of the order parameter. Fol- 
lowing the convention of Ref. E| the scaling function 
is raised to a negative exponent, with <f> = and 
(3 = 0.326. The solid curve is the result obtained in 
this work for the WR mixture; the dashed curve is the 
HCSW result of Ref. fej Also shown is the exact small x 
limiting form y = 1 + x/2. 

FIG. 2: Order parameter of the WR mixture as func- 
tion of the distance from the critical point. The dashed 
line has a slope (3 — 0.326, corresponding to the 3D Ising 
exponent. 

FIG. 3: Scaling function of the coexistence diameter 
for the WR mixture. Open circles show simulation re- 
sults obtained using the FSS algorithm of Ref. [H The 
dashed curve is a fit to the simulation data using the 
approximant of Eq. J5J ■ Also shown is the exact small x 
limiting form y — x/2. 

FIG. 4: Coexistence diameter of the WR mixture as 
function of the distance from the critical point. Open 
circles are simulation results obtained in this work using 
the FSS algorithm of Ref. ^3- The dashed curve is a 
fit to Eq.lQ]). The black dot marks the critical density 
obtained from the fit, where the vertical line indicates the 
uncertainty. The arrows mark estimates of p cr reported 
in Ref. H (JEA) and Ref. (Buhot), where the vertical 
lines again indicate the uncertainty. The inset shows k 
as function of t. Open circles are simulation results; the 
dashed curve, which essentially overlaps the simulation 
data, is a three-parameter fit of the form n — a\t^ a + a-2, 
with fit parameters a, and a. 

FIG. 5: Cumulant analysis of the WR mixture. Shown 
is the first order cumulant U± as function of the fugacity 
zb for various system sizes L as indicated. The inset 
shows the slope of the cumulant Y± at the critical point as 
function of L. All data were obtained along the symmetry 
locus za — z-q- 
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